Auxiliary Field quantum Monte Carlo for Strongly Paired Fermions 
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We demonstrate that the inclusion of a BCS importance function dramatically increases the 
efficiency of the auxiliary field method for strong pairing. We calculate the ground-state energy 
of an unpolarized fermi gas at unitarity with up to 66 particles and lattices of up to 27'^ sites. 
The method has no fermion sign problem, and an accurate result is obtained for the universal 
parameter ^. Several different forms of the kinetic energy adjusted to the unitary limit but with 
different effective ranges extrapolate to the same continuum limit within error bars. The finite 
effective range results for different interactions are consistent with a linear term proportional to 
the Fermi momentum times the effective range. The new method described herein will have many 
applications in superfluid cold atom systems and in both electronic and nuclear structures when 
pairing is important. 



The study of strongly interacting Fermi systems is one 
of the central themes and major challenges in physics. 
Superfluidity in unpolarized cold atomic Fermi gases, 
which has been demonstrated both experimentally and 
theoretically, provides a prototypical example. The ex- 
perimental ability to use a Feshbach resonance to adjust 
the strength of the potential between the atoms allows 
an exploration of the physics over many length scales. A 
particularly interesting regime is at unitarity where the 
scattering length diverges and the range of the potential 
can be neglected. Since the particle density provides the 
only length scale, the ground-state energy _Eo is propor- 
tional to the free fermi gas energy Epci 
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The ability to quantitatively understand the proper- 
ties of this system represents a great triumph of many- 
body physics. Many experiments and calculations have 
been performed for the unitary Fermi gas. Initial qual- 
itative agreement was found between theory[TJ [2] and 
experiment[3]-[S]. More precise recent experiments have 
yielded C = 0.39{2)[B] and 0.41(1)[7], with smaller val- 
ues obtained very recently by Zwierlein, et al.[5]. Fixed- 
node Diffusion Monte Carlo (DMC) calculations [D 0111- 
[T2] have always included a Bardeen-Cooper-Schrieffer[T3] 
(BCS) trial wave function to guide the Monte Carlo walk 
and provide the fixed node constraint |14] needed to over- 
come the fermion sign problem. As is well known, these 
calculations provide an upper bound, with the current 
best value ^ = 0.383(1) [III ES]- 

In this paper we show that exact calculations can 
be performed to accurately determine the ground-state 
properties of the unpolarized Fermi gas. A new method is 
introduced to allow the use of a BCS trial wave function 
in the auxiliary-field quantum Monte Carlo (AFQMC) 
approaches of Zhang and coworkers [El [16]. Using the 
new approach, we perform calculations with several forms 
of the kinetic energy term that all give the correct con- 



tinuum limit but with different finite effective ranges to 
study the convergence with particle number and lattice 
sizes, and to obtain the dependence of ^ on the effective 
range. An exact result is obtained for the value of £,. 

Quantum Monte Carlo simulations play a key role in 
addressing the challenge of strongly interacting Fermi 
systems. The AFQMC method has been applied to a 
variety of systems in several fields. With equal numbers 
and masses of up- and down-spin fermions and an attrac- 
tive interaction, there is no fermion sign problem. The 
formalism presented here allows the use of a BCS impor- 
tance function, which drastically improves the efficiency 
in this situation. In general applications, a sign or phase 
problem is present, which is controlled by a constraint, 
also using the importance function [T^l [TB] . Hartree-Fock 
or free Fermi gas (FG) type of importance functions have 
typically been used. This approach has been shown to be 
very accurate in many condensed matter models and op- 
tical lattices ^7\ , quantum chemistry [THl , and solid state 
materials [THj. Now BCS importance functions (or anti- 
symmetrized geminal power (AGP) in chemistry) will sig- 
nificantly improve our ability to deal with the sign prob- 
lem in systems where pairing is important, and enhance 
the capabilities for quantum simulations in strongly cor- 
related systems in general. 

The AFQMC method, in both-zero[2n] and finite- 
temperature formulations [21], has also been applied to 
the unitary Fermi gas. Precise results require simultane- 
ously large lattices, so the system is dilute, and a large 
number of particles for an accurate approach to the con- 
tinuum and thermodynamic limits. Many such calcula- 
tions have been Derformed |20l [22ll25j . but the variance of 
the method limited the results to relatively small number 
of particles and lattice sizes so that, as we demonstrate 
below, the results are unlikely to have converged to the 
thermodynamic limit. 

The range of the van der Waals interaction in cold 
atoms is small (e.g. about 3 nm in ^Li[26 ) compared 



to interparticle spacing so that the short range structure 
of the interaction is unimportant; the results are com- 
pletely determined by the form of the kinetic energy and 
the scattering length. For an N^ lattice, the equivalent 
Hamiltonian is 
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Here ipjs is the destruction operator for a fermion of spin 
s on lattice site at position Tj. For odd lattice sizes 



with 



2N, 
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operators are Cfe 



1, the k are given by "^{rij-x + n^y + n^i) 
< n^.n^i^n^, < Nc- The k space destruction 

— 3/2 Y~> ^—ik-ri 

k 



Ly, , 



Nr'^Ei 



'ipjs, and the den- 



sity operators are riis = ipl^ipis- 

To reach the continuum limit, we need to take the limit 
of zero particle density, p = N/N^ — )■ 0, in the context of 
the Hubbard model (i.e., replace L by N/^). Equivalently, 
because of scale invariance, we can think of the system as 
a discretized representation of a supercell with fixed size 
L, and take the fc-space cutoff to infinity. In cither case, 
we then take the number of particles, N, to infinity. In 
this limit only the behavior of e^ for fc <C — is important, 
where a = L/Nf^ is the lattice spacing. Thus a variety 
of kinetic energy forms can be used as long as they are 
quadratic for k much smaller than the cutoff. In this 
work we present results for 
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The superscript 2 and 4 indicate the highest power of fc, 
while e'-''-' is the Hubbard model hopping kinetic energy 
offset by a constant so that it is zero at fc = 0. 

For two particles, the Hamiltonian is separable, and 
the solution of the Lippmann-Schwinger equation for low- 
energy s-wave scattering gives the phase-shift equation, 
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where V indicates the principal parts integration, and the 
k space sums are cut off by the lattice spacing a. The 
effective range expansion is 



k cot (5o = —a 



(5) 



where a is the scattering length and r^ the effective range. 
Since we are interested in the unitary limit, we adjust 
U to have a~^ = 0. Both t\. and e^ ' have nonzero 
effective ranges. The extra parameter /3 in ej^ can be 
adjusted to make the effective range zero. The values for 
the parameters are given in table ll] 



Energy 


TJ 2ma.-' 


13 


r.a-' 


4"' 


-7.91355 


- 


-0.30572 


4^' 


-10.28871 


- 


0.33687 


ei^' 


-8.66605 


0.16137 


0.00000 



TABLE I. The parameters that give infinite scattering length 
for two particles in an infinite lattice for the various kinetic 
energies. The /3 value for the ej. kinetic energy has been 
adjusted to give zero effective range, r^. 



The AFQMC algorithm uses branching random walks 
to project the ground state from an initial trial state with 
the imaginary-time operator exp[— iJr]. Because the in- 
teraction is attractive, there is no fermion sign problem 
for equal numbers of up and down fermions studied here, 
and no path constraint is required. A walker is a set of 
N single-particle orbitals. Initially, the orbitals for the 
up spin particles are taken to be identical to those for 
the down spin particles. The two-body interaction term 
is broken up using a Hubbard-Stratonovich (HS) trans- 
formation which has only positive weights, and treats the 
up and down spin particles identically. Therefore, the up 
spin orbitals remain identical to the down spin orbitals 
during the propagation. We will show below that the 
usual form for a singlet paired BCS trial function also 
gives no fermion sign problem. 

The walker states are given by specifying the orbital 
coefficients. These can be specified on the real space lat- 
tice (j)nj or as momentum space coefficients (j)n,k related 
to each other by a discrete Fourier transform. If we begin 
with real orbitals on the real space lattice, the orbitals 
remain real when propagated. The momentum space or- 
bitals therefore satisfy (pn.-k = 4>n k- ^^^ orbitals are 
orthonormalized at each step. This is needed to limit 
roundoff error, but the mathematical expressions are cor- 
rect without it. The orthonormal orbitals therefore sat- 
i'^fy Z]fc ^,* ,fc<^ni,fc = ^nm, and the corresponding opera- 
tors, Wns = Z]fe<^n,fcCfcs, Satisfy {Wns,wl^^,} = SnmSss'- 

The walker state is 

N/2 



m = l[ ^It^ljo) 
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Because iV <^ N^ and because the imaginary-time his- 
tory of the walk need not be retained, this formalism is 
much more efficient than the usual lattice formulations 
for the ground state of dilute gases. 

Using a discrete HS transformation[27 , the potential 
energy propagator is 

^ " W}=±i 
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where tanh u — ~ tanh (^^) ■ The kinetic energy prop- 
agator is 



Gri^t) =exp 



k 



At 
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Given a choice of one of the N^ set of fields, the ap- 
pUcation of the Trotter breakup of one term of the prop- 
agator on a walker \W) gives another walker \W') times 
a weight w'{{(j}, W) that depends on the set of HS vari- 
ables {a} and \W), 

G{{<j},At)\W) ^ GT{f)Gv{{a},At)GT{f)\W) 

^w'i{a},W)\W'). (9) 

The propagation above consists of (1) Multiply each 
4'n,k by exp(— iffcAt). (2) Fast Fourier transform to 
obtain 0„.j in real space. (3) Multiply each 4>n.i by 
exp (2ucri — ^C/At). (4) Fast Fourier transform to ob- 
tain 4)n,k in momentum space. (5) Multiply each (j)n.k by 
exp(— 2efcAt). (6) Update the weight from non-operator 
terms. 

Including importance sampling reduces the fluctua- 
tions, by changing the sampling so that it is non-uniform, 
without biasing the results. We want to sample walkers 
\W) from {'^T\W){W\il;{t)) where 
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The contribution of a walker \W) at the next time step 
is then 
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We want to sample the set of HS variables {<j}, from 
the unnormalized probability distribution given by the 
square brackets. The normalization which, to order At^ 

is the local energy expression g-ts ['^^(^)+^^(^')l-^r)^*, 
will give the weight of the sampled walkers. Once we 
have sampled the {a} values, the new normalized walker 
is given by the second term of Eq. 11 We make sure 



that regions where the trial function is small are sampled 
adequately to eliminate trial-function bias. 
The particle projected BCS state is 
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where fk = Vk/uk in the usual notation. The overlap of 
a walker with the BCS state is 



{W\BCS) = 
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The contraction needed is Wns'c\g = 4'n,k5ss'- The two 
creation operators in the BCS pair must be contracted 
with some 'Wm-\ and some Wn^, giving a term 



Anm = WnlWm^ ^ /fe4t'^ 



t it 



-ki 



k 



= 2_^ (t'n-kfk<t>m,k — 2_^ 4>n,kfk4>m,k ■ 
k k 

Taking all the different possible full contractions, 
(Wipes') =det^. 



(14) 



(15) 
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matrix A are the An 



where the elements of the 
of Eq. [14| 

The overlap, Eq. |15[ is positive when, as in the stan- 
dard singlet paired BCS solutions, fk > 0. We write 
det^ — det [SB'''] where B^ is the hermitian conjugate 



matrix of B and the matrix elements of the 



N 



Nl 



2 '^ ^^k 

trix B are _B„fc = (l)n,k^/~fk ■ Applying the Cauchy-Binet 
theorem, each of the determinants of the submatrices of 
B is multiplied by the determinant of the corresponding 
hermitian conjugate submatrix. The determinant of A is 
a sum of positive terms, so that our BCS trial function 
gives no sign problem. 

Two estimates of the energy are used, the growth en- 
ergy just measures the growth of the weights in the ran- 
dom walk. The local energy can be calculated using con- 
tractions similar to those above. Other observables can 
be calculated similarly. A complete derivation for a gen- 
eral BCS state will be published elsewhere. Here we give 
the result 



El{W) 



{W\H\BCS) 



2tr [A-^C] 



{W\BCS) 
+ UY,{ir[A-'E{q)\ir[A-^E\q)\ 

-XT{A-^E[,q)A-^E\q)\ + tr [^-^^^(q)] } , 
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where Aim(q) = J2k^n.k+qfk4>Tn,k+q, Cn„i = 
J2k4>n,k'^kfk<t>ni,k, and Enrniq) = T,k K,k+qfk<f>rn,k- 

The matrix elements of D and E are convolutions which 
are efficiently computed with fast Fourier transforms. 
The computational cost of using the BCS I'^t) is sim- 
ilar to using a single Slater determinant. 

We have calculated the ground-state energy for dif- 
ferent particle numbers and lattice sizes. The time-step 
errors have been extrapolated to zero within statistical 
errors, and walker population biases have been checked 
and were found to be negligible for the population sizes 
used. The imaginary time step is « 0.01-0.05 Ep^ , the 
total propagation time is of order 10-30 Ep^ and 2,000- 
20,000 random walkers are used in the simulations. 

For A^ = 4, we found that the use of BCS importance 
functions reduced the statistical error by a factor of 10, or 
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FIG. 1. (color online) The calculated ground state energy 
shown as the value of ^ versus the lattice size for various 
particle numbers and Hamiltonians. 



for a variety of lattice sizes and find agreement with Ref. 
[25] . The error bands in the figure provide least-squares 
estimates for the one sigma error based upon quadratic 
fits to the finite-size effects. The fits are of the form 
E/EpG = C + ^P^^^ + Bp^/^. For the interactions tuned 
to Te = 0, a fit with A fixed to zero is used. Including 
a linear coefficient in the fit yields a value statistically 
consistent with zero. 

The extrapolation in lattice size for the fc^ and Hub- 
bard dispersions show opposite slope as expected from 
the opposite signs of their effective ranges. The extrap- 
olation to p — >■ is consistent with ^ = 0.372(0.005) in 
all cases. Our final error contains statistical component 
and the errors associated with finite population sizes and 
finite time-step errors. This value is below previous ex- 
periments, but more compatible with recent experimental 
results of the Zwierlein group |8J. 



100 X reduction in computer time, compared to the usual 
FG importance function. The improvement increased to 
1500 X for iV = 38 in a 12'^ lattice. For larger systems, 
the discrepancy is much larger still; indeed the statistical 
fluctuations from the latter are such that often meaning- 
ful results cannot be obtained with the run configurations 
described above. 

In Fig. [T] we summarize our calculations of the energy 
as a function of p^^^ where p = N/N^, and the particle 
number is N = 38, 48 or 66. We plot ^, Eq. [l] where we 
have in all cases used the infinite system free-gas energy 



EpG^iV^ with k% = 37r2^ 



as the reference. 
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Hamiltonian 


iV 


^ err 


A err 


4^' 


14 
38 
66 


0.39 0.01 
0.370 0.005 
0.374 0.005 


0.21 0.12 
0.14 0.04 
0.11 0.04 


4^' 


38 
48 
66 


0.372 0.002 
0.372 0.003 
0.372 0.003 




4'^' 


4 
38 
48 
66 


0.280 0.004 
0.380 0.005 
0.367 0.005 
0.375 0.005 


-0.28 0.05 
-0.17 0.03 
-0.05 0.03 
-0.13 0.03 



TABLE II. Energy extrapolations to infinite volume, zero 
range limit for various particle numbers A^ and different 
Hamiltonians. The term linear in the effective range, A, is 
also shown where it is not tuned to zero. 

DMC calculations have found converged results when 
using 66 particles [TTl [T^ . and our results confirm this. 
The differences between 38 and 66 particles are rather 
small. Our calculations with 14 particles show a signif- 
icant size dependence, and with 26 particles the effects 
are still noticeable. These are not shown on the figure. 
We have also computed the energy for 4 particle systems 



FIG. 2. (color online) The ground-state energy as a function 
of kpr^: comparison of DMC and AFQMC results. Dashed 
lines are DMC results, shifted down by 0.02 to enable com- 
parison of the slopes. 

We have also examined the behavior of the energy 
as a function of kpre for finite effective ranges. It has 
been conjectured [28j that the slope of f is universal: 
C(''e) = S, + Skprg. Of course a finite range purely attrac- 
tive interaction is subject to collapse for a many-particle 
system, but a small repulsive many-body interaction or 
the lattice, where double occupancy of a single species is 
not allowed, is enough to stabilize the system. Our re- 
sults are consistent with the universality conjecture. In 
particular our results for zero effective range approach 
the continuum limit with a slope consistent with zero. 

Figure b^ compares the AFQMC results for the e\. ' in- 
teraction with the DMC results JTj V2\ for various values 
of the effective range. The AFQMC produces somewhat 
lower energies than the DMC, consistent with the upper- 
bound nature of the DMC calculations. For the slope S of 
^ with respect to finite r^, the fit to the TV = 66 AFQMC 
results yields S = O.I I (.03). Similar fits to the AFQMC 
data with the Hubbard dispersion el, for iV = 66 yield 



a linear term of 5* = 0.12( 
with the DMC results of S 



.03). Both are in agreement 
= 0.12(.01).[29] 



In summary, we have shown how to incorporate a 
pairing importance function into auxiliary field quantum 
Monte Carlo algorithms and used it to treat the uni- 
tary Fermi Gas. This algorithm, for attractive interac- 
tions and equal spin populations, is exact and can be ex- 
tended to large lattices and strong interactions. We find 
^ = 0.372(.005) using a variety of interactions tuned to 
unitarity. We also find a slope of the ground state energy 
with effective range of S= 0.12 (.03) for the different lat- 
tice and continuum Hamiltonians. This method should 
be useful without modification for the entire BCS/BEC 
transition and for studying many other properties of cold 
Fermi gases. It can also be applied to a wide variety of 
problems in other strongly-correlated fermions, in areas 
ranging from cold atoms to condensed matter to quan- 
tum chemistry to nuclear physics. 
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